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Abstract 

We study low-speed flows of a highly compressible, single-phase fluid in the pres- 
ence of gravity, for example in a regime appropriate for modeling recent space-shuttle 
experiments on fluids near the liquid-vapor critical point. In the equations of motion, 
we include forces due to capillary stresses that arise from a contribution made by strong 
density gradients to the free energy. We derive formally simplified sets of equations in 
a low-speed limit analogous to the zero Mach number limit in combustion theory. 

When viscosity is neglected and gravity is weak, the simplified system includes: 
a hyperbolic equation for velocity, a parabolic equation for temperature, an elliptic 
equation related to volume expansion, an integro-differential equation for mean pres- 
sure, and an algebraic equation (the equation of state). Solutions are determined by 
initial values for the mean pressure, the temperature field, and the divergence-free part 
of the velocity field. To model multidimensional flows with strong gravity, we offer 
an alternative to the anelastic approximation, one which admits stratified fluids in 
thermodynamic equilibrium, as well as gravity waves but not acoustic waves. 
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1 Introduction 



Near the liquid- vapor critical point, many of the thermophysical properties of a fluid exhibit 
a singular behavior. For instance, the isothermal compressibility and the isobaric thermal 
expansion coefficients, as well as the isobaric specific heat, all diverge strongly at the critical 
point. Critical enhancement effects are also encountered in the behavior of the thermal 
conductivity and the viscosity in the vicinity of the critical point, while the thermal diffusivity 
approaches zero. These singularities play a major role in the thermal equilibration of near- 
critical fluids. 

Understanding the effect of singular fluid properties on dynamics is not always straight- 
forward. For example, it has been shown that even though thermal diffusivity is small, 
temperature changes in a near-critical fluid can occur rapidly via a mechanism which causes 
adiabatic pressure changes in the bulk of the fluid fZ5| . This adiabatic mechanism cre- 
ates a strong coupling between temperature changes occurring at the fluid boundaries and 
the temperature response in the interior of the fluid. It works as follows. A temperature 
perturbation applied at the boundary of a fluid causes an expansion in the fluid near the 
boundary. Through the medium of sound waves, this produces an adiabatic pressure change 
in the interior of the fluid, and a consequent change in the temperature, much more rapidly 
than could be accomplished by thermal diffusion acting alone. Near the critical point, where 
the fluid becomes highly expandable and compressible, the adiabatic mechanism dominates 
the early thermal response and creates a 'critical speeding-up' phenomenon. This critical 
speeding-up has been observed in earth-bound and low-gravity experiments ||, |5| . 

In contrast to the short time-scale of the thermal response, experimenters have observed 
a much longer time-scale for the equilibration of density variations [BO, IT5| . Near the critical 



point, the divergence of the compressibility and the influence of gravity can create strong 
macroscopic density gradients (upon which microscopic density fluctuations are superim- 
posed). Although early adiabatic processes act rapidly (within seconds) to accomplish most 
of the temperature changes, most of the relaxation of the density distribution to a new equi- 
librium state is a non-adiabatic process driven by the much slower (hours-long) process of 
heat diffusion. 

Recently, Boukari, Pego, and Gammon |J studied the combined effects of the adiabatic 
mechanism and earth's gravity on the equilibration process in near-critical xenon, using a 
system of equations which includes not only the adiabatic effect, but also one- dimensional 
fluid motion and heat advection. In numerical simulations of a temperature step exper- 
iment, they found that the onset of a transient diffusive regime occurs about ten times 
sooner than estimated by Onuki, Hao and Ferrell [25| in the zero-gravity case, due to the 
generation of a small temperature gradient by the adiabatic pressure quench in the pres- 
ence of gravity. Boukari et al. also observed that over periods of many hours, no single, 
exponentially-decaying mode was ever observed to dominate the diffusive equilibration pro- 
cess. This conclusion is consistent with measurements and computations of Zhong and Meyer 
and Kogan, Zhong and Meyer ||17|| . 

The conclusions drawn in these works were based on results derived from one-dimensional 
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systems of equations that do not account for multi-dimensional flows. It is not yet clear how 
multi-dimensional flows affect equilibration under gravity. Zappoli et al. [35|] have performed 
computations of two-dimensional buoyant flow for a van der Waals fluid fairly near the 
critical point (IK above critical) where stratification effects are not yet very large, using an 
'acoustic filtering' technique resembling the approach taken in this paper, and have observed 
an unusual convection pattern. In the zero-gravity situation, R. F. Berg || has pointed 
out certain differences between one-dimensional and corresponding three-dimensional results 
regarding the late diffusive regime. 

The purpose of this paper is to systematically derive systems of multi-dimensional equa- 
tions which govern the dynamics of a low-speed, highly compressible, single-phase fluid in 
the presence of gravity. We shall include forces due to capillary stresses that arise from a 
contribution made to the energy by strong density gradients. Although there are no sharp 
interfaces between phases in equilibrium in the one-phase regime just above the critical tem- 



perature, strong gradients can be generated as transients [[L4]]. We anticipate that nonlocal 
effects related to the energy of 'diffuse interfaces' could play a big role in generating transient 
flows in certain circumstances that are accessible to experimental observation. 

Our starting point is the general hydrodynamic equations expressing the conservation 
of mass, linear momentum and energy for a compressible fluid with heat conduction and 
gravity. We suppose that the equation of state is appropriate for conditions near the critical 
point, and presume that the fluid is in local equilibrium. This implies that the time-scale 
of interest is longer than a local relaxation time and that the critical point is not so near 
that the correlation length is macroscopic. These assumptions appear to be reasonable for 
describing the regimes studied in recent experiments. 

To account for the influence of steep density gradients on energy, we adopt a constitutive 
structure compatible with that described in the work of J. E. Dunn and J. Serrin 0. Dunn 
and Serrin permit the constitutive quantities (such as the Cauchy stress tensor) to depend 
upon the gradient of the density as well as upon density and temperature, so that capillary 
effects can now be included in the equations. We account for such effects in the simplest 
way, adding a squared gradient term to the Helmholtz free energy density. The resulting 
system of equations in section 2 can be used to study the influence of hydrodynamics on 
heat transfer on a time scale appropriate for resolving sound waves. 

Here, however, we are interested in relaxation and flow phenomena that occur very slowly 
compared to the time it takes for sound waves to cross the spatial domain. A systematic pro- 
cedure for obtaining simplified equations for compressible flows on long time-scales was intro- 
duced by Rehm and Baum [pjj and by Majda and Sethian [^TJ in the context of combustion 
theory. In section 3 we apply this procedure to the system at hand. We non-dimensionalize 
the equations and estimate the size of the various terms, taking parameters appropriate to 
a typical experiment of interest. To obtain a simplified set of equations, we neglect terms 
which make the smallest contribution compared to the other terms. 

When the effect of gravity is weak, the result is a simplified model in which the leading- 
order pressure is spatially constant, and which accounts for multi-dimensional fluid motions 
influenced by the effects of thermal expansion and contraction, gravitational compression, 
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thermal diffusion, viscosity and capillarity. With viscosity included, this simplified model 
consists of equations (|3.12| )- (|3.16| ) below, plus the equation of state. 

When viscosity is neglected, the simplified system includes a hyperbolic equation (for 
velocity), a parabolic equation (for temperature), an elliptic equation (related to volume 
expansion), an ordinary integro-differential equation (for mean pressure), and an algebraic 
equation (the equation of state). In a forthcoming work ||, we prove that the simplified 
model equations are evolutionary; i.e., we show that solutions are determined by suitable 
initial data. The initial data required to determine the solution consist of the temperature 
field, the mean pressure, and the divergence-free part of the velocity field. The density is 
determined from the temperature and mean pressure via the equation of state. 

Earth's gravity creates a strongly nonlinear density profile in equilibrium when temper- 
ature is close to the critical point, due to the high compressibility. This is not well-modelled 
by the system ( |3.12|) -( |3.16|) . In section |4] we study how to model multidimensional flows 
with strong gravity. If we enforce hydrostatic balance at leading order, as was done for the 
one-dimensional case in ||, we find that generally the flow must remain strictly stratified, 
and any vertical fluid motions are due only to horizontally uniform density changes. 

In order to admit nontrivial convective flows or gravity waves, one can assume that 
entropy is constant at leading order and neglect heat conduction. This corresponds to the 
assumptions made by Ogura and Phillips [23 in deriving 'anelastic' equations for atmospheric 
circulations. Here we obtain anelastic equations valid for a general equation of state. The 
assumption of approximately constant entropy may not be compatible with thermodynamic 
equilibrium near the critical point, but this assumption may be relevant for describing some 
experiments that have been performed at near-constant entropy in order to reduce the effect 
of density gradients @, p2| . Because heat conduction is neglected, though, the anelastic 
equations do not contain the fast adiabatic heat-transfer mechanism mentioned above. 

We find a better alternative if we return to the weak-gravity scaling and make a simple 
modification of the momentum equation. We retain the effect of the pressure correction 
on density in the gravity force term, in the spirit of the Boussinesq approximation. This 
modification does not change the formal validity of the approximation. But the new system, 
consisting of equations ( |4.18|) - ([4.20|) below, correctly captures strongly stratified thermo- 
dynamic equilibria, admits multidimensional flow including gravity waves but not acoustic 
waves, and includes heat conduction and the adiabatic mechanism. 



2 Basic hydrodynamic equations 

To start, we consider the general hydrodynamic equations expressing the conservation of 
mass, linear momentum, and energy for a compressible fluid with heat conduction and grav- 
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ity: 

% - -,V-v, (2.1) 

p Wt = T:Vv - Vc i- ( 2 - 3 ) 

Here T is the Cauchy stress tensor, p is the density, v is the velocity, e is the specific 
internal energy density, — g is the gravitational acceleration, q is the heat flux, and D/Dt = 
d/dt + v • V. 

To account for the influence of density gradients on energy in a manner compatible with 
the second law of thermodynamics, we apply a theory for Korteweg-type fluids described by 
J. E. Dunn and J. Serrin |J. A recent review of related theories and their applications to 
diffuse-interface modeling has been given by Anderson, McFadden and Wheeler |T|. Dunn 
and Serrin replace the energy balance equation with 

De , , 

p— = T : Vv — V • q + V ■ u (2.4) 

J J L 

where V • u represents a contribution made to the power by capillary effects due to the strong 
density gradients. They then postulate that e, T, q, u, as well as if), the specific Helmholtz 
free energy density, and s, the specific entropy density, are given by constitutive relations 
that depend only on the pointwise values of p, T, Vp, VVp, VT, and Vv, where T is the 
temperature. 

In order to guarantee that the equations of motion are compatible with the second law 
of thermodynamics as expressed by the Clausius-Duhem inequality, Dunn and Serrin deduce 
that the constitutive relations must be such that 

^(p,T,|V/f), s = -j£, e = ^-T^, (2.5) 

and 

u= -p(V • v)mVp + (2.6) 

where m = 2p(dip/dM)(p, T, M) with M = |Vp| 2 . Furthermore, in the case without viscos- 
ity, the stress must take the form 

-P 2 tt- + PV • (raVp) j 1 - mVp® Vp. (2.7) 
dp J 

In (|2.6|) , the quantity u measures the "static" part of the capillary work flux u. For a class 
of materials including those we shall consider, V • u = 0, so uo has no effect on the energy 
balance equation and can be ignored. 
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We choose the simplest possible form for the specific Helmholtz free energy density ip 
that is consistent with Dunn and Serrin's theory. Namely, we suppose that m is constant, 
and that 

A rji 

^ = ^T) + -\Vp\ 2 . (2.8) 
2p 

This yields an expression for total free energy that appears in van der Waals' theory of 
capillarity [p8 ], for example. From (|2.5| ) it follows that 

e = e(p,T) + |-|Vp| 2 , s = -^(p,T), (2.9) 

where e = ip — T(dip/dT). We assume that the heat flux is given by q = — kVT, with 
thermal conductivity coefficient n = k(p,T). And we suppose that the stress is given by 
( ^.7| ) with the addition of Newtonian viscosity terms, so that 

T= (^-p+ — \Vp\ 2 + mpApj 1 -mVp® Vp + A(V- v)l + 2pD, (2.10) 

where the pressure p is given by 

p = p(p,T)=p 2 ^(p,T), (2.11) 

D = |(Vv + Vv T ) is the symmetric part of Vv, and the viscosity coefficients have the form 
A = A(p, T), p = p(p, T). For compatibility with the Clausius-Duhem inequality ||, one 
requires that A + |p > 0. 

When equations ( |2.6|) , ( |2.10|) , (|2.8|) and (|2.9| ) are substituted into the balance laws (|2.1|) , 
(f2.2[), (|2.4|), the result is the following system of equations for a viscous, heat-conducting 
fluid with capillary stresses: 

P— + Vp = -pg + mpVAp + V(AV ■ v) + V- (2/xD), (2.12) 



Dt 



DT ( \ c v \ Kt Dp , ,_, 



Dt \ c p J a p Dt 



+ (pcp) _1 V ■ (kVT) 



H-(pCp)- 1 (2pD : D + A(V- v) 2 ) , (2.13) 



Here J^t = p 1 (dp/dp)r is the isothermal compressibility, a p = — p 1 (dp/dT) p is the isobaric 
thermal expansion coefficient, and 

d£ dp 
c v = c v (p, T) = — , c p = c p (p, T) = c v + p~ x a v T— 

are the specific heat capacities at constant volume and at constant pressure, respectively. 
The system ( f^.l2j )-( ]2.14| ) differs from the standard system for a viscous, heat-conducting 
fluid by the addition of the single term mpVAp in the momentum equation. 
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3 Equations for slow flows with weak gravity 



3.1 Scaling for the flow regime 

In order to see how the system of equations ( |2.12|) — ( p.!4j ) can be appropriately simplified, we 
non-dimensionalize and scale the variables in a manner appropriate to a typical experiment, 
as follows. The critical density, critical pressure, and critical temperature are denoted p c , 
p c , and T c , respectively. Let x a be a characteristic length, t a a characteristic time, and 
v a = x a /t a . Let c Pa , c Va be characteristic specific heat capacities at constant pressure and 
at constant volume respectively, and let T = c Pa /c Va . Let K a be a characteristic thermal 
conductivity, and let X a , p a be characteristic viscosity coefficients (for convenience we assume 
A a = Pa)- Let g a = x a /t 2 a , and introduce the following non-dimensional variables and 
constants (non-dimensional quantities are indicated by an asterisk, * ): 

P = PcP\ p = p c (l+p aP *), T = T c (l + T a T*), 
x = x a x*, t = t a t\ v = v a \\ g = g a g*, 

C v = c v a c *,i c p = c p a C *pi K = K a K *) ^ = A a A*, p = Pap* i 

t p * \dp* ) Tt ' p p* \dT* 

Here, p a and T a are non-dimensional scales characterizing the deviation of the pressure and 
temperature from critical. One tries to choose the scales so that the dimensionless variables 
p*, T*, v* and their derivatives are of the order of 1 in the flow regime of interest. 

The non-dimensional equations corresponding to equations ( |2.12| )-( |2.14| ) for a viscous 
fluid are 

p*^- = -M~ 2 Vp* + cp*VAp* - p*g* 



Dt* 



+Re _1 (V(A*V ■ v*) + V ■ (2//D*)) , (3.1) 



Dt* \ c*J a* Dt* p*c* 



S[ ^-B* : F>* + —(X7- v*) 2 ] , (3.2) 
\p*c* p*c* p K ' I 



Dv* DT* 

K*-f- = al— — -V-v*. (3.3) 
1 Dt* p Dt* 



where the gradient and the divergence are taken with respect to the non-dimensional spatial 
variable, and D/Dt* = d/dt* + v* • V. 

The dimensionless constants M 2 , Dt, c, Re -1 , and S are defined by 

M 2 = P_A n =J^, c =!^, Re- 1 = ^4, S ^ 



PcPa Pc c p a X "a X a Pc X a Pc c p a TcT a ta 
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The parameter M is proportional to the Mach number v a /c s , where the sound speed is 
given by c 2 = (dp/dp) s = c p /c v pK T . Re is the Reynolds number, D T is a non-dimensional 
diffusivity, and c is a non-dimensional coefficient of capillarity. 

A typical flow regime in which we are interested is one considered in the paper by Boukari, 
Pego, and Gammon |J. The fluid is xenon, with critical parameters 

p c = 1.11 x 10 3 ^|, Vc = 5.84 x 10 6 Pa, T c = 289.72 K. 
nr 

A typical experimental cell is about half a centimeter in radius, so we choose x a = 10~ 3 
m. The time scale of interest ranges from a fraction of a second to hours. For now we take 
t a = 1 s and postpone further discussion. The effect of earth's gravity tends to become 
important within about 30 mK of T c , so we take T a = 10~ 4 . In this temperature range 
and near the critical density, an appropriate model equation of state is the restricted cubic 



model, with coefficients for xenon (see |22[ and the references therein) . From this model, as 
in we find it appropriate to take 

Pa = 6T a , c Pa « 3.3 x 10 6 — I— , T- 1 « 1.8 x 1(T 4 . 

kg - K 

We have not introduced a separate scale for deviations of the density from critical because 
these can be rather large, of the order of 10%. 

We estimate a characteristic value for thermal conductivity as in ||, using the approach 



described in [31, 29 1 to calculate the background term and divergent part near the critical 



point. To estimate the viscosity, we use the results of |3l], Table III], also see |pq| . As a 
result we find it appropriate to take 

«; a «2.5xl0 _1 — -, p a « 5 x 10~ 5 Pa ■ s. 

s • m ■ K 

Finally, we estimate the constant capillarity coefficient m using the power law represen- 
tation m ~ /xt derived by Rowlinson and Widom [28], where the correlation length £ and 



susceptibility xt = p 2 Kx at the critical density are given |30| by the power laws 



£~£o|ATT", Xt ~ ^^\AT*n. (3.4) 

Pc 

Here AT* = (T — T c )/T c = T a T* w 10~ 4 , and we use from |3(| the critical exponents 7 = 1.19 
and v = 0.63, and for xenon in the range T > T c approximately £0 = 1-9 x 10~ 10 , C = .0813. 
Using these values we obtain the estimate 

(5.84 x 10 6 )(1.9 x 10" 10 ) 2 , in 4x 07 A in 18 

m w ^— ^— -^(10" 4 )" 07 «4x 10 -18 . 

(1.11 x 10 3 ) 2 (.0813) V ; 

From the estimates above, we obtain the following estimates for the dimensionless pa- 
rameters: 

M 2 ^3xl0- 7 , D T ^6xl0^ 5 , c^4.4xl0~ 3 , 
Re" 1 w 4.5 x 10" 2 , S ^5x 10~ 13 , 
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and we find 



l3 <? , /3p" 



9.81x10', -|~1, ^] ~1. 



The sound speed c s ~ 80 m/s. 

Regarding these parameters, several points are worthy of comment. First, note the effect 
of considering longer time scales. As t a increases, S decreases and M 2 decreases quadratically, 
Dt and Re -1 increase, and c and g* increase quadratically. Second, experiments performed 
in low earth orbit are reported to experience typical accelerations of 10~ 4 to 10~ 6 times 
earth's gravity 22| ; this would make g* of order 1. (Another way to obtain g* of order 1 is 



to consider a faster time scale like t a = .01 s.) Also, we note that the value of c becomes 
of order 1 when the characteristic length is replaced by a capillary length x C8iP for which 
1 = t 2 a mp c / 'a; 4 ap . For our flow regime, we estimate x cap ~ x a (4.4 x lO -3 ) 1 / 4 ~ 260 microns. 
In the recent ZENO experiment, observations were performed using light scattering through 
a fluid layer 100 microns thick JT3 . 



We should also comment on the effect of proximity to the critical temperature. As T a 
approaches zero, we have seen that the capillary coefficient diverges very weakly, with expo- 
nent 7 — 2v —0.07. The viscosity and the sound speed also diverge at a very slow rate, 
changing only modestly over the experimental range of interest. The specific heat c v also 
diverges weakly, with exponent —a ~ —0.11. The compressibility, thermal expansivity and 
specific heat c p = c v + xtT p~ 3 (dp / dT) 2 diverge strongly, all with exponent —7 ~ —1.19. 
The thermal conductivity diverges less strongly, like c p /£/i |JD], p. 22], with approximately 
the exponent —(7 — v) « —0.56. So we see that as T a approaches zero, none of the nondi- 
mensional constants above has a very strong dependence on T a , though we can expect the 
nondimensional diffusivity Dt and S to decrease. 



3.2 Reduced equations for slow flows 

To obtain simplified equations for compressible flows, we proceed formally in a manner 
motivated by the treatments of Rehm and Baum [[27] and Majda and Sethian ||21|| . Since we 



are interested in longer time scales, we regard M and S as small and let po, T , v denote 
assumed asymptotic limits of p*, T*, v* respectively, as M 2 and S are taken to zero. In this 
process we regard g* as fixed and of order one, corresponding to a low-gravity environment 
for our flow regime. Multiplying the momentum equation ( |3.1| ) by M 2 and taking M 2 to 0, 
we require 

Vp (xV*) =0, (3.5) 

and therefore po = po(t*) is constant in space. From (|3~2| ) -(|3T3l ) , the asymptotic equations 
for temperature and pressure are 

DT ( „_ x c v q \ K to dp D T 

1-T — — + V-(k VT ), (3.6) 



Dt* \ CpoJ a p0 dt* p Q c p0 

dp DT 

'dr - ap0 m*~ 



dp DT 

K to-j- = u p0 — -V-v . (3.7) 
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Here the asymptotic density p (x*,£*) = p*(p (t*),T Q (x.*,t*)) from the equation of state, 
and elsewhere the subscript indicates a nondimensional coefficient that is evaluated at 
(p ,T ), e.g., K (x*,t*) = K*(p ,T ). Equation ( |3.7| ) implies the mass continuity equation 
Dp /Dt* + p V • v = 0. 

Next we derive a reduced momentum equation. We suppose the flow occurs in a bounded 
domain Q with v = on the boundary dQ. We use the fact that every square-integrable 
vector field v has a unique orthogonal decomposition of the form 

v = w + V0, where V ■ w = and w ■ n\ dn = (3.8) 

(n is the outward unit normal to dQ). We write Pv = w, so P denotes the orthogonal 
projection of square-integrable vector fields onto solenoidal vector fields. Observe that if 
Pf = 0, then f = V/i for some function h, and conversely. Applying P to the momentum 
equation ( |3.1| ), then, eliminates the term of order M~ 2 . Taking M 2 to produces 

~ Po T^ ~ PoS * + cpoVApo + ^ ( V ' (MVv + Vv£)) + V(A V • v )) 



= P 



where p = n*{po,T Q ), A = A*(p ,T ). 

From this we infer that there must be a scalar function pi = pi(x*,t*) such that the 
expression above in brackets equals Vpi, that is, 

Po^ + Vpi = -p g* + cpoVApo + ^ (V ■ (p (Vv + Vv^)) + V(A V ■ v )) . (3.9) 

Note that at this point we do not presume that the pressure p* ~ p + M 2 pi to order M 2 . 
Klainerman and Majda [16[ have found that at order M 2 there is an acoustic correction to 



pressure that depends on fast time and space scales. 

Equations (pi.6|) , (|3.7|), and ( ^.9| ) are the simplified equations for a viscous fluid. In order 
that our approximations be self-consistent, we must require that at the initial time, as M 2 
and S tend to zero we have 

p*(x*, 0) -> p (0), T*(x*, 0) ^ T (x*, 0), v*(x*, 0) ^ v (x*, 0). 
3.3 Reformulation 



The full system ( |2.12| )-( |2.14| ) is appropriate for describing compressible fluid flow on acoustic 



time scales. In the flow regime for xenon that we have described, such time scales are short, 
since the sound speed is of the order of tens of meters per second. The simplified system 



( p.6|) , (|3.7p , (|3.9|) represents an 'acoustic filtering' of the full system that describes flow on 
time scales that are long compared to acoustic. The pressure is maintained spatially constant 
through a process mediated by sound wave propagation. (An asymptotic description of this 
process for a near-critical van der Waals fluid was given by Zappoli and Carles ]34| in one 
dimension with no viscosity.) 
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As is easy to check, both the full system and the simplified system are compatible with 
the Clausius-Duhem inequality 

Ds „ /q\ 

We next reformulate the simplified system (|3.6| ) , ( p.7|) , (|3.9|) into an equivalent form which 
better reveals its evolutionary character. This reformulation will be necessary for our future 
purpose of analyzing the initial-value problem. We will omit the subscripts and superscripts 
and return to dimensional quantities for notational convenience, writing 

P = Pcpo, P = p c (l + PaPo), T = T c (l + T a T ), v = v a v . 

We employ the decomposition v = w + V0 described in ( |3.8|) . Substituting the temper- 
ature equation ( PTB| ) into the pressure equation ( |3.7| ) and solving for V • v, we obtain 

c v Kt dp a n , mX 

V-v = A0 = — —-7- + — V • (kVT). (3.10) 

This elliptic equation has a solution with V0 • n = v • n = on dfl (and then V0 is uniquely 
determined), if and only if a solvability condition holds, namely 

c v Kt dp i a p 



+ -^V • (kVT) )dx = 0. (3.11) 



Solving this equation for dp/dt, we get 

dP rr M In ( a p/P c p)^ ■ ( kVT ) dx /« 10 X 

— = H[t).= — . (3.12) 

at J n {c v K T /c p ) ax 

This integro-differential equation is the pressure evolution equation. The remaining equa- 
tions of the system can be written as 

DT f l-^)^H(t) + —V-(KVT), (3.13) 



Dt \ c p J a p pc p 

p— = - Vtt - pg + cpVAp - P^f 1 

+V • (//(V(w + V0) + V(w + V0) T )) + V(AA0), (3.14) 

A0 = -2^1 H (t) + ^V-(KVT), (3.15) 

V-w = 0, (3.16) 

along with the equation of state p = p(p, T). (Here 7r = p c p a M 2 pi.) 

From ( |3.15| ) evaluated at time t = 0, we obtain the compatibility condition 



A0| 



t=o 



Cp PCp 



(3.17) 



t=o 
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This equation imposes a constraint on the gradient part V0(x, 0) of the initial velocity field 
v(x, 0), and leaves freedom for the choice of the initial solenoidal component w(x, 0), as long 
as V • w(x,0) = 0. 

As we show in || , to solve the initial- value problem for these equations, it is appropriate 
to specify initial data for the leading-order pressure, temperature, and divergence-free part 
of the velocity field: 

p(0) = Po, T(x, 0) = T (x), w(x, 0) = w (x) (3.18) 
where V • wo = 0. Initial data for the velocity will take the form 

v(x, 0) = v (x) = w (x) + V0 o (x), (3.19) 
where V</>o is determined by the compatibility condition ( |3.17| ). 



3.4 The adiabatic time scale 



We conclude this section by indicating how the simplified equations ( |3.12|) -( |3.16|) , in nondi 



mensional form, can be used to very roughly estimate the time scale ti for the adiabatic 
mechanism described in the introduction to produce a rapid bulk temperature response to 
boundary heating. 

The nondimensional form of ( |3.12| ) is 



dp = TD T J n (a P o/poc p0 )V ■ (K VT ))dx 
dt* f n (c v0 K T0 /c p0 ) dx 

We consider a homogeneous fluid initially at equilibrium, whose boundary temperature is 
raised rapidly. As the boundary temperature is changed, a thin thermal boundary layer is 
created next to the wall. The width of this layer increases with time through heat diffusion, so 
is roughly given by y/t*Dr using (|3.6|). In the boundary layer, we may roughly approximate 



the nondimensional temperature Tq by a function of the form f(s/\^t*Dr), where s is the 



distance to the boundary. Treating the coefficients in (|3.20|) as constant, we estimate AT[ 



o 



(t*Dx) -1 f"(s/yj t*D-r) in the boundary layer. The integrand in the numerator is then of order 
(t*Dx)~ l in the boundary layer and zero elsewhere, and the integrand in the denominator 
is of order 1. Suppose the fluid domain is a cube with dimensional side length L = x a L*. 
Then the order of dpo / dt* is given by 



dp QL^VFD^^Dt)- 1 { 3QT 2 D t \ 1/2 

dF ~ VDt = [ ~L*r 1 ■ (3 - 21) 



The time integral of this expression produces an order one change in p (hence in Tq), when 
t — t\ — t*t a = L 2 T~ 2 /144D, where D = n a /p c c Pa is the characteristic thermal diffusivity. 
For the fluid parameters corresponding to the flow regime which we have described above, 
in a cell with side length L = 10~ 2 m we estimate t\ ~ 3.3 x 10~ 4 s. 
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This time is longer but not much longer than the acoustic time t c for a sound wave to 
cross the cell, given by t c = L/c s pa 1.25 x 10 -4 s. It is unrealistic to expect, however, that 
the boundary temperature can be raised so rapidly in experiment. So what this estimate 
indicates is that for processes in which the boundary temperature varies slowly compared 
to the acoustic time, the boundary-to-bulk coupling provided by the integro-differential 
equation for the mean pressure is efficient in effecting bulk temperature changes. 

The above estimate for t\ is consistent with the results of Onuki and Ferrell [J24|, except 



for the geometric factor of 144 appearing in the denominator. We expect that diffusion 
dominates the equilibration at approximately the time = L 2 /lAAD pa 10 4 s. Onuki, Hao 
and Ferrell |25[ characterize the intermediate regime between the long times t/t\ > T 2 and 
the short times t/t\ = 0(1) by the geometrical mean t- mt /ti = T so that t int = Tti. For our 
flow regime we estimate ti nt ~ 2 s. 



4 Multi-dimensional flows with strong gravity 
4.1 Motivation 

In this section, we re-examine the equations of motion in the case of strong gravity. Mo- 
tivating us is the problem of describing near-critical fluid flows and equilibration in earth's 
gravity. Recall that |g*| ~ 10 4 in the flow regime considered in section |3.1| with earth's grav- 
ity. The simplified system fl3.12|) -( [3.16|) fails to capture some key features of equilibration 



in this situation. 

In equilibrium, temperature is constant and density is stratified according to the basic 
equation of hydrostatic balance, 

Wp = -pg, (4.1) 

and the equation of state. (We will neglect capillarity in most of this section.) Denoting 
equilibrium temperature by T e and density by p e (z), the equilibrium density gradient satisfies 

^r(z) = -XT(p e (z),T e )g. 
az 

As temperature approaches T c , the critical temperature, XT{.Pc,T e ), the susceptibility on the 
critical isochore, diverges as in Q3.4| ). Thus the density gradient develops a singularity at the 
level of critical density, and the density profile becomes highly nonlinear. In Fig. 1 we plot 
density profiles for xenon in equilibrium at 1G, using the restricted cubic equation of state 
as in 0. 

For the simplified system ( |3.12| )— (|3.16f) , however, the density in equilibrium is constant, 
given by p = p(p e ,T e ). The hydrostatic balance W = — pg from ( |3.14| ) can be interpreted 
as supplying an 0(M 2 ) correction to the leading order pressure. This correction is linear in 
z and could be used to generate a density correction (by linearizing the equation of state, 
for example). But as it stands, the system ( |3.12|) -( |3.16| ) relies only on the leading-order 
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x 10~ 3 

Figure 1: Density vs height for near-critical xenon at 1G. 
T - T c in niK: 1 (dash), 10 (solid), 100 (dash-dot), 1000 (dot) 

density to determine the fluid flow. One can therefore expect this system to be inaccurate in 
modeling phenomena such as deep convection and gravity waves in near-equilibrium states. 

In order to model related phenomena, researchers studying small-scale atmospheric circu- 
lations frequently approximate the continuity equation with the 'anelastic' continuity equa- 
tion 

V • (pv) = 0, (4.2) 

where p(z) is usually defined either as the density in an adiabatically stratified, horizon- 
tally uniform reference state, or as the horizontally-averaged actual density. Batchelor || 
introduced an equation equivalent to (|4.2|) . The name 'anelastic' was given by Ogura and 
Phillips | Jffi| , who derived ( |4.2| ), together with approximate momentum and thermodynamic 
equations, through a systematic scale analysis. Important assumptions in their analysis are 
that: (i) all deviations 59 of the 'potential temperature' from some constant mean value 9 a 
are small (this is equivalent to a similar statement for entropy variations), and (ii) the time 
scale of the disturbance is similar to the time scale for gravity wave oscillations. The terms 
neglected in their approximation are formally an order e = 59 /9 a smaller than those which 
are retained. Their anelastic system does not support sound waves, does support gravity 
waves and conserves energy. Ogura and Phillips define p as the density in an adiabatically 
stratified, horizontally uniform reference state. 

Some of the largest errors in the Ogura-Phillips anelastic approximation are reported 
to be generated by large deviations of the mean state potential temperature (or entropy) 
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from a constant reference value. Several authors (Dutton and Fichtl [ll|, Wilhelmson and 



Ogura p2| , Lipps and Hemler [[0|], Durran |T(|) have presented alternative 'sound-proof 



equations in which thermodynamic variables associated with the adiabatic reference state 
are replaced with horizontally-uniform averaged or approximate actual values. These authors 
make different approximations in the momentum equations, but (except for Durran) they all 
obtain a continuity equation of the form ( }4.2|) in which p is defined as a horizontally-averaged 
approximation to the actual density. 

In the remainder of this section, we describe three possible models for multi-dimensional 
flows under strong gravity, making different assumptions about how to balance and approx- 
imate terms in the nondimensionalized system (|3.1|)-(|3~3"D: 



(1) At first we scale so as to enforce hydrostatic balance at leading order. We find then 
that to be consistent, at leading order the thermodynamic variables must remain hor- 
izontally uniform. Also the vertical velocity must remain horizontally uniform, unless 
we assume the entropy is spatially constant to leading order. 

(2) If we indeed assume that entropy is spatially constant to leading order and also neglect 
heat conduction, we obtain a generalization of the anelastic approximation valid for a 
general equation of state. 

(3) Neither choice so far admits both thermodynamic equilibrium states and non-stratified 
multidimensional flows. We find, however, that if we return to the scaling adopted 



in section |3.1| , and modify the momentum equation so that the gravity force term 
incorporates a pressure correction self-consistently, then we get a system of the same 
formal accuracy that admits physically correct equilibrium states exactly. Moreover 
the new system supports gravity waves but not acoustic waves. 

4.2 Strongly stratified flow 

Because equilibria are governed by the equation of hydrostatic balance, if gravity is strong 
it is natural to try and balance the pressure gradient term with the gravity term in the 
nondimensional momentum equation ( |3.1|) . Thus we regard M~ 2 and |g*| to be of the same 
order. This same approach was taken in [H for one-dimensional flows. 

In this approximation, we rescale the gravitational acceleration, writing g* = M _2 g 
where g = 0(1). It will prove instructive to replace (p,T) by (p, s) as the thermodynamic 
state variables, where s is the specific entropy density. In terms of these variables, we write 
the equation of state as p = p(p, s). In nondimensional form, we write p* = p*{p*, s*), where 
entropy is nondimensionalized via the relation s = s a s*. It is convenient to take s a = c pa T a ; 
this will be discussed further in the next subsection. 

Starting from the nondimensionalized system ( |3.1| )-( |3~3| ) with the replacement g* = 



M g, we proceed as in section and regard M and S as small, obtaining the temperature 



evolution equation ( |3.6| ) and continuity equation ( |3.7| ). In terms of the leading order pressure 
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p and entropy s , these equations are equivalent to 



Po (l + T a T )^ = D t V-(k VT ), (4.3) 
Dt* 

1 Dp a p0 ( Ds \ , v 

p (l + T a T )— — - V ■ v . (4.4) 



p c z s0 Dt* poc p o V -D** 

The leading-order terms in the momentum equation now just yield the equation of hy- 
drostatic balance, 

Vp = -Pog- (4.5) 

This equation imposes tight restrictions on the spatial dependence of the leading-order ther- 
modynamic state variables. In particular, it is necessary that po is a function only of height 
z* and time t* and is independent of the horizontal variables (x*,y*). Then from ( |4.5| ) it 
follows that p — Po{z*,t*) is also horizontally uniform. The equation of state now implies 
that so = so(z*,t*) as well. The entropy and pressure equations (fO)-(fO) now read 



1 fdpo dp \ a p0 d f dT \ 

w ottt = Dt— [ k -V-v . (4.7) 



PoCso \ 9t* 9z* J PoCpo dz* \ dz* 

Here Wq is the vertical component of the velocity v . 

If we assume dso/dz* ^ 0, then to satisfy ( [Op consistently, the vertical velocity must 
be horizontally uniform: Wq = Wo(z*,t*). We can eliminate the horizontal components of 
velocity from ( |4.7| ) by integrating over (x*,y*); if the cell walls are vertical or periodic we 
can express the result in terms of density as the usual one- dimensional continuity equation, 

£ + ^ = * (4 . 8) 



The system consisting of the three equations (f4.5|), ( |4.6| ), ( |4.8| ) is equivalent to the one- 



dimensional system (14)-(16) in || for temperature and pressure, which was expressed using 
a Lagrangian variable z' = J* po(h, t*) dh in place of the height z* £ [0, L\. With t' = £*, in 
present notation this system takes the form 

p (z',t') = p a {t')-gz\ (4.9) 
dT ( c v0 \ K T0 dp a D T d ( dT \ 

dp a f™Po 1 a po (dT /dtf)dz' 

dt' " f^Po'KTodz' ■ l4 - iiJ 
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(Here m = J Q L po(h,t*) dh is the total linear mass.) This system determines the evolution of 
the temperature, pressure and vertical velocity. Vertical fluid motions in this model are due 
solely to density changes that occur in a horizontally stratified manner — the equations do 
not support acoustic or gravity waves. 

The horizontal components of the momentum equation at order 1 govern the horizontal 
fluid flows within each material layer corresponding to z' =const. Presuming that p\ yields 
an 0(M 2 ) correction to the leading order pressure, we may write the horizontal components 
of the momentum equation as follows: Let u denote the horizontal components of velocity 
v and the subscript h denote differentiation with respect to the horizontal variables = 
(x*,y*). Then the horizontal momentum equation for u(x/i, z' ,t') is 

Po (j£ + (u ■ V fc )uJ + V h p x = ^ (^ A h u + po^y (a^o|^ J ) • (4.12) 



Horizontal motions are incompressible, because • u = from ( [4.7|) and (O). When 



viscosity can be neglected, then the equations corresponding to different fluid layers decouple, 
and fluid layers can exhibit arbitrary independent two-dimensional incompressible flows. 

The vertical component of the momentum equation at order 1 can be used in deter- 
mining higher-order corrections of order M 2 for horizontally averaged density, pressure and 
temperature. For the sake of brevity we omit further discussion. 

4.3 The anelastic approximation 

In order to admit convective flows with nontrivial vertical circulation and/or gravity waves, 
we can suppose that the entropy is spatially constant at leading order. As we have indicated, 
this is the same as the assumption made by Ogura and Phillips that potential temperature 
variations are small compared to a reference value. 

At leading order, then, the fluid density, pressure and temperature are in hydrostatic 
balance and adiabatically stratified, with a (negative) adiabatic temperature gradient. We 
must neglect heat conduction to maintain entropy constant; this restricts the time scale, and 
any flows generated will be adiabatic. 



To describe flows, we use the notation of section |3J] for nondimensional and leading- 
order quantities, and nondimensionalize entropy according to s = s a s* where s a = c pa T a 
(as discussed below). We postulate that the nondimensional pressure p* and entropy s* are 
given to order 0(M 2 ) by 

p* ~ Po (z*) + MV(x*, t *), s*~s + M 2 Sl (x*, O, 

where so is a constant. (This neglects any fast-time acoustic corrections that may be present 
as discussed in [0.) Then we expect the density p* ~ Pa{z*) + M 2 pi(x.*, £*), where from the 
equation of state, 

* = -^ + (l + T.T )*!SV (4.13) 
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Here the nondimensional coefficients are evaluated at (po,s ). 
The momentum equation at order 1 is 



Since heat conduction is neglected, the entropy s% is convected with the flow, satisfying 

^ = 0. (4.15) 
Dt* y ' 

Because dp /dt* = 0, the leading-order continuity equation yields 

V • (p v ) = 0. (4.16) 



This constraint on the velocity should be used in solving (|4.14 ) to determine p±. Note 



however, that p± will not be completely determined by the constraint — Given any solution 
of ( |4.14j ), pi can be replaced by adding an arbitrary time-dependent multiple of a solution 



to the linearized hydrostatic balance equation 



9 ~ 



dz* c 2 s0 

The equations ( [4.14| )-( |4.16| ) correspond to the anelastic equations of Ogura and Phillips, 
generalized for an arbitrary equation of state. 

As a model for slow flows of fluids near the critical point, these equations have some 
drawbacks: First, heat conduction is entirely neglected, so the fast adiabatic mechanism for 
rapid thermal response is not accounted for within this model. Also the effect of fluid flow 
on thermal relaxation cannot be evaluated. Moreover, describing states in thermodynamic 
equilibrium is problematic. At rest, the equations permit an arbitrary horizontally uniform 
entropy correction s\. It may or may not be consistent with the derivation of the equations 
to take this to correspond to the equilibrium entropy profile (meaning the hydrostatic profile 
at constant temperature). 

Regarding this point, we can ask, what is the size of the nondimensional equilibrium 
entropy gradient in the regime of interest? To estimate this, we need to identify a typical 
entropy change in a process of interest. Consider a fluid at equilibrium at one temperature, 
subject to a temperature change at the boundary of order T c T a . During the early development 
of the thermal boundary layer, we may suppose roughly that the entropy change in the 
boundary layer occurs at constant pressure, so it is of order s a = c pa T a since (ds/dT) p = c p /T 
and T a is small. Then it seems reasonable to nondimensionalize the entropy by letting 
s = s a s*. 

Now, the equilibrium entropy gradient satisfies 

ds I ' ds\ dp c p — c v ( dT\ 
Tz = \dp~J T Tz = ~^T~ \ dp) p 
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Nondimensionalizing this expression, we find that up to a quantity of order one, 



^^_M 2 |g*| = -^^. (4.17) 

dz* p c p a 

Under the specific conditions considered in section |37I] , ds* /dz* ~ —.003. This is moderately 
small but not quite of order M 2 ; it would make ds\/dz* of the order of — 10 4 . Also, since 



p a ~ 6T a , the last member of ( f4.17|) shows that ds* /dz* diverges as temperature approaches 
critical, like T~ l as T a — > 0. 

It seems that in our flow regime, the entropy gradient may not be small enough for 
the anelastic approximation to be reasonable for states near true equilibrium (at least at 
the level of critical density). But states near constant entropy do possess experimental 
interest. Physically, the condition that an inviscid stratified fluid is stable against convection 
is ds/dz < |P-8|| . In order to obtain 'pseudo-equilibrium' states with near- uniform density 
profiles, ground-based experiments have been suggested on near-critical fluids in which a 
steady-state heat flux is maintained so as to achieve a marginally stable entropy profile |[22|| , 
also see ||. Note that the adiabatic density gradient satisfies p~ l dp/dz = —g/c 2 , which 
is about 1/500 m _1 in our flow regime, so the density is close to constant in centimeter- 
sized cells. Provided that heat conduction can be neglected on the time scale of interest, 
the anelastic approximation could be useful to describe flows near such 'pseudo-equilibrium' 
states. 



4.4 Slow flows including equilibrium 

In trying to enforce hydrostatic balance at leading order, we have found that either flows 
remain strictly stratified, or equilibrium states are not admissible as entropy must be constant 
at leading order in M 2 . In this section we return to the scaling as it was done in section 
f3.2| , in which the leading-order pressure turns out to be spatially constant. Working in the 
spirit of the Boussinesq approximation, we can decide to selectively retain some terms of 
higher order in M 2 where it would be most useful, without affecting the formal accuracy of 
the system. 

We propose no alteration in the equations ( |3.6| )-( |3T7| ) for leading-order temperature and 
mean pressure. (But now po and T will depend on M 2 , through the coupling to the mo- 
mentum equation.) To order M 2 , supposing that the pressure is given by p* ~ po + M 2 pi, 
the approximation to the density p ~ po = P*(po,2~o) can be improved to 

p^p = p*(p + M 2 Pl ,T ). 

We propose to make this improvement only in equation ( |3.9| ) for the velocity. The new 
system of equations governing leading-order temperature, pressure and velocity (neglecting 
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capillarity) is 



if = ( 1 - r "r)^ + 7r v ' (K ° VT ° ) - (418) 

Lit \ Cpoj a p o at PoCpo 

dp DT 

K ™^f = «po^-V-v , (4.19) 

P^ + ^Pi = -pg* + ^(V-(/io(Vvo + Vv^)) + V(AoV-v )). (4.20) 

The temperature evolution and continuity equations (|4.18|) -( |4.19| ) may be replaced by the 
equivalent pair (|4.3| )- (|4.4| ), in which the material derivative Dpo/Dt* can be replaced by the 
ordinary derivative dpo/dt*. 

As with the anelastic equations, the pressure correction pi should be determined in 
solving ( 4.20|) to satisfy the implied constraint on the divergence of velocity from ( 4.19| ). 



The pressure correction will not be unique, but any two solutions p\ and p\ that correspond 
to the same (p ,T , v ) will be related by 

Vpi + p*(p + M 2 Pl , T )g* = Vpx + p*(p + M%, T )g*. (4.21) 

Therefore the difference p\ — pi is a function of z* and t* that is determined solely in terms 
of a function of t* by solving an ordinary differential equation. 

The system consisting of equations fl4.18|) -( |4.2C| ) can be reformulated to better reveal 
its evolutionary character exactly as in section |3.3| . In dimensional form with pressure 
p(t) + n(x,t), where p = p c {l + p a Po) and n = p c p a M 2 pi, one obtains exactly the system 
([Tg)-(gT§) except that ( gTg ) is replaced by 

+V • (/i(V(w + V0) + V(w + V0) T )) + V(AA0), (4.22) 

where the (now dimensional) density p is determined from pressure and temperature by the 
equation of state: p = p(p + 7r,T). We anticipate that, like system ( |3.12j )-( ]3.16| ), solutions of 



( |4.18| )- (|4.20| ) are determined by initial values for the temperature field, mean pressure, and 



divergence-free part of the velocity field. 

The system ( |4.18| )- (|4.20| ) admits as rest states true equilibrium states with constant 
temperature and hydrostatic balance between the total pressure p + M 2 pi and the improved 
density p = p*(p + M 2 pi,T Q ). It filters acoustic waves but admits gravity waves, as we shall 
show in the subsection to follow. Moreover, heat conduction need not be neglected, so the 
adiabatic effect can be modeled. It should be interesting to study what flows are generated 
when large density changes in the thermal boundary layer are present together with very 
stable equilibrium entropy profiles away from the boundary. 

We remark that it is evidently not necessary to replace po by p in the acceleration 
term of ( |3.9| ) to gain true equilibria as rest states. In some circumstances it may be more 
convenient not to make this replacement. But it turns out that equation ( |4.20| ) is slightly 
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more convenient when we study the linearized equations for gravity waves; see the next 
section. 

The system ( |4.18| )-( |4.20| ) certainly has shortcomings. First, in the flow regime described 
in section |g*| = gt 2 a /x a is still rather large. This problem diminishes if a larger space 
scale or smaller time scale is relevant. The system can be expected to lose formal validity if 
solutions become large or singular, as may well happen in a nonadiabatic convection process. 
Another point is that while the leading-order total energy is conserved in time for the weak- 
gravity equations (|3.12|) - (|3.16|) , this is not strictly true for the system ( |4.18| )-( [4.20|) . The 
time derivative of the leading-order total energy is formally of order M 2 instead. 



4.5 Gravity waves in the linear approximation 

We wish to verify that in the linear approximation near a stably stratified rest state for 
which entropy decreases with height, the system ( 4.18| )- (|4.20|) admits gravity waves but not 



acoustic waves when heat conduction and viscosity are neglected. We shall also show that, in 
the special case of a perfect gas with an exponential density profile and constant sound speed, 
in the limit of large wave number the gravity-wave frequency approaches the Brunt- Vaisala 
frequency iV corresponding to a compressible fluid. This frequency satisfies 

N 2 = -1**L - ? (4.23) 

Pe dz Cj 

The second term does not appear in the usual treatment of gravity waves for a stratified 
incompressible fluid, in which the density is advected with the flow, cf. |33| . 

When heat conduction and viscosity are neglected in the system ( |4.18[ )-( ^.2U| ), the leading- 



order pressure is constant in time as well as space. Consequently the temperature (and en- 
tropy) are advected with the flow, and the velocity field has zero divergence. In dimensional 
form with p = p + tt and p = p(p, s), the governing equations take the form 

w t - °- < 424 > 

V-v = 0, (4.25) 

Dv 

P^ + Vp = -pg. (4.26) 

Near a rest state where (p, s, v) = (p e (z), s e (z), 0), we write 

P = Pe+P, s = s e + s, v=(u,v,w) 
and p e {z) = p{p e , s e ). Then we linearize, obtaining 

ds . , N 

— + ws e {z) = 0, 



V-v = 
dt ^ \\dp J r ' \ds 



Pe^r + Vp = -[IT) P+(^r) S 



p 
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where the coefficients are evaluated at (p e ,s e ). Let us suppose that periodic boundary 
conditions are specified in the horizontal variables (x, y) . We look for normal modes with 

(p,s,u,v,w) = e ^ kix+k2y -^ (p, s,u,v, w) (z) , 

and eliminate p, s and the horizontal components of velocity from the system. Note that 
the sound speed satisfies c~ 2 = (dp/dp) s , and define 

1 dp e , . If dp\ ds e , , g 



a{z) = -j^ l3[z) = -j\-£) ii =a[z) -fr (4 - 27) 

(The last identity holds due to hydrostatic balance.) Then for the vertical velocity component 
w(z) we obtain the equation 

d 2 w dw / 9 k 2 Bq\ 

+ P-T+[k 2 qJi)w = 0, (4.28) 



dz 2 dz \ a; 2 

where k 2 = k\ + k\. The vertical velocity must vanish at the top and bottom boundaries. 
So, given a horizontal wave number k, possible oscillation frequencies u are determined by 
solving the eigenvalue problem in ( |4.28| ) with Dirichlet boundary conditions. 

We note that in the special case of perfect gas at constant temperature, the density profile 
is exponential and sound speed is constant, and so a(z) and (3(z) are constant. Then ( 4.28 ) 



has explicit solutions of the form w(z) = e' 3z l 2 sin(nz), whence the gravity-wave dispersion 
relation is given by 



k Bq 

- 2 = < 4 - 29 > 



In the limit k 2 — ► oo, this expression approaches (3g = N 2 , where N from (|4.23|) is the 
Brunt-Vaisala frequency for a compressible fluid. For comparison, for a fully compressible 
fluid in which one starts with Dp/Dt + pV ■ v = in place of ( [4.25| ), the dispersion relation 
in this special case is 



«» = | (e + n 2 + y ± y (* 2 + » 2 + t) 2 - 4fc2 f 

For large k 2 , the plus sign yields uj 2 ~ c 2 k 2 , corresponding to acoustic waves, and the minus 
sign yields uj 2 ~ iV 2 , corresponding to gravity waves. The magnitude of any oscillation fre- 
quency uj that satisfies (|4.29f) is less than N, showing that the system (|4.24| )- (|4.26| ) supports 
gravity waves but not acoustic waves. 

In general, when the coefficients in ( |4.28| ) vary with z, we can obtain an upper bound on 
oscillation frequencies as follows. Let q(z) = exp 

(- J* P(0 rf C), multiply equation ( ggg ) by 
qw and integrate over z £ [0, L], from bottom to top. One obtains 

k 2 r L 

q(z)w'(z) 2 dz+— (uj 2 - (3g)q(z)w(z) 2 dz = 0. 

(jJ .In 
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Since the second integrand cannot be everywhere nonnegative, it follows that 

uj 2 < max(3(z)g. (4.30) 



The dispersion relation in ( |4.29| ) is very similar to that obtained in the usual case of an 



incompressible fluid when one assumes the density is advected with the flow. Starting from 
the governing equations 

f t = 0, ,4.3!) 

Vv = 0, (4.32) 

P^ + Vp = -pg, (4.33) 

one finds in similar fashion that the equation corresponding to ( |4. 28| ) is 

d 2 w dw f 2 k 2 ag 



+ oi-j- + [k — )w = 0, (4.34) 



dz 2 dz \ uj 2 
and for an exponential density profile the dispersion relation is 



U = k 2 + n 2 + \a 2 - (435) 



In the limit k 2 — > oo this approaches ag = Nq, where N is the usual Brunt- Vaisala frequency 
for an incompressible fluid. In many circumstances the difference between N and N may be 
negligible, but it is interesting that the dispersion relation arising from the system ( f4.24p - 
( |4.26| ) more faithfully approximates the compressible case in this respect. 



4.6 Final remarks 

We close with a few remarks intended to clarify the differences between the new system 
( |4.18| )- (|4.20| ) and the anelastic system ( |4.13| )-( [4.16| ). For purposes of comparison, we neglect 
heat conduction. In this case, the new system (|4.3| ), Q4.4|) , ( [4.20|) becomes 



^ = 0, (4.36) 

Dt* v 1 

V-v = 0, (4.37) 

^ + VP1 = -^* + ^( V -(/ i o( Vv o + Vv o)) + V ( A oV-v )), (4.38) 

where p = p*(po + M 2 p 1: s ). 

For each system, the source of the constraint on velocity (( f4.37p or ( [4.16[ ) respectively) 
is the continuity equation 

^ + p V • v = 0, (4.39) 
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where p = P*(Po> s o)- F° r the anelastic system, the leading order entropy s Q is constant 



and the pressure p is a function of z determined by the hydrostatic balance equation (|4.5|) , 
so Dp /Dt* = Wodpo/dz* . Using hydrostatic balance we can also write ( |4.16| ) in the form 
essentially given by Batchelor 0]: 

V • v - w -^- = 0. (4.40) 

For the new system (}4.36|) -( }4.38|) , since p is a constant when heat conduction is neglected, 
we have Dp /Dt* = and this is why the velocity field is divergence- free. Note that if we 
seek to 'improve' equation ( }4.39|) by replacing p by p = p*(po + M 2 pi, s ), then we recover 
the original fully compressible system without simplification! Indeed, the essential difference 
between the fully compressible system ( |3.1|) - (|3.3|) (neglecting heat conduction and viscous 
power terms), and system fl4.36|) - ([4.38|) with p* = po + M 2 pi, is precisely that in the new 



system a term proportional to M 2 Dpi/Dt* is neglected in the continuity equation. 

This point suggests a modification to the system (|4.36| )- fl4.38|) in a situation with possible 



relevance for atmospheric circulations. Suppose gravity is rather strong but the pressure 
correction pi does not happen to deviate significantly (more than 0(1)) from some time- 
independent reference state p(z*) that determines a reference density profile p(z*) via a 
hydrostatic balance equation 

Vp = -p(z*)g?. (4.41) 

Then we replace po in ( |4.39| ) by po = p*(po + M 2 p, s )- Note that po can depend on (x*, y*, t*) 
as well as z* through so- Since M 2 Dp/Dt* = M 2 wo(dp/dz*), the continuity equation be- 
comes 

V • v - wo J--?- = 0, (4.42) 

C s0 PO 

where the coefficient c 2 is evaluated at (p + M 2 p, s ) and we have used M 2 g* = g. This 
equation replaces ( [4.37| ), without changing the formal validity of the approximation. 
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